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Abstract 

A molecular-level model is used to study the mechanical response of empty cowpea chlorotic 
mottle virus (CCMV) and cowpea mosaic virus (CPMV) capsids. The model is based on the 
native structure of the proteins that consitute the capsids and is described in terms of the 
atoms. Nanoindentation by a large tip is modeled as compression between parallel plates. Plots of 
the compressive force versus plate separation for CCMV are qualitatively consistent with continuum 
models and experiments, showing an elastic region followed by an irreversible drop in force. The 
mechanical response of CPMV has not been studied, but the molecular model predicts an order of 
magnitude higher stiffness and a much shorter elastic region than for CCMV. These large changes 
result from small structural changes that increase the number of bonds by only 30% and would be 
difficult to capture in continuum models. Direct comparison of local deformations in continuum 
and molecular models of CCMV shows that the molecular model undergoes a gradual symmetry 
breaking rotation and accommodates more strain near the walls than the continuum model. The 
irreversible drop in force at small separations is associated with rupturing nearly all of the bonds 
between capsid proteins in the molecular model while a buckling transition is observed in continuum 
models. 
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INTRODUCTION 



Understanding the properties of self-assembled molecular structures such as virus cap- 
sids is of great current interest. Capsids are conglomerates of proteins that protect virus's 
genomes. They are generally strong mechanically, especially when they envelope tightly 
packed DNA at high internal pressures. For instance, in bacteriophage 029, the pressure is 
of order 50 atm [ll-y]. Capsids that protect single stranded RNA have an order of magni- 
tude lower pressures, as evidenced by Michel et al. [Si] for the cowpea chlorotic mottle virus 
(CCMV), but they are still quite sturdy. 

One way of assessing the elastic properties of capsids is by making nanoindentation 
measurements with an atomic force microscope. In particular, it has been found 5| that 
the CCMV capsids with their RNA removed resist nanoindentation less (by about 30%) 
than their full counterparts. Both types of CCMV capsids are highly elastic - the initial 
Hookean reversible regime persists up to nanoindentions varying between 20% and 30% of 
the diameter and apphed forces of 600 pN and 1000 pN for the empty and full capsids 
respectively . These forces are noticeably larger than characteristic forces needed for 
a mechanical unraveling of proteins. As an example, titin, which is a reasonably strong 
protein, requires a force of about 204 pN 7, ^ to unravel whereas the particularly strong 
scaffoldin protein c7A requires a force of about 480 pN j^. 

The sturdiness of the capsids comes from both the mechanical properties of the individual 
proteins and the way the proteins are bound together. Capsids have a variety of shapes and 
sizes, but most have highly symmetric icosahedral structures. Early models treated the 
capsid as a homogeneous spherical shell, first in the thin shell limit {4, 5| and then including 
the finite width of the shell and nonlinear response fll|. The latter work showed that a 
Hookean model for finite width shells did not reproduce the linear force displacement curve 
observed in experiments, while standard nonlinear models for rubberlike materials did give 
a linear response. The icosahedral symmetry of viruses has been included in the thin shell 
imit ay , la, llO| and using the mass density to determine a spatially varying shell thickness 
1^ . but with homogeneous elastic response. 

Many of the above models are roughly consistent with experimental behavior. The force 
rises linearly with displacement for an extended range and then drops suddenly, signaling 
an instability. The linear region can be fit by assuming elastic moduli with a reasonable 
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magnitude, but the fact that different models can achieve this suggests that the result is 
not sensitive to variations in local moduli and density that may be present in capsids. 
For example, one may expect that the number and geometry of bonds within proteins 
are different than those between proteins, and this heterogeneity is difficult to capture in 
continuum theories. Another issue is that the instability at higher forces in continuum 
models is associated with a buckling instability without any breaking of molecular bonds. 
The latter may be important in real capsids, but is difficult to parameterize in continuum 
models. More detailed models are needed to assess potential variations in elastic properties 
within the shell and the potential role of bond breaking. 

All-atom simulations are an appealing alternative, but are limited to relatively short 
indentation times of order Ins compared to experimental times of 10ms or more. Studies of 
mechanical properties have included an all-atom analysis of the low frequency vibrations in 
the satellite tobacco necrosis virus |l3l and an elastic network model analysis of eigenmodes 
in CCMV [l^ and HK97 viruses [l5, 16|. Recent simulations by Zink and GrubmuUer 
have examined indentation of the southern bean mosaic capsid using an atomistic treatment 
with surrounding water. The force on a small (~ Inm) sphere representing the tip was 
measured as a function of the distance from the center of mass. They find a rather short 
elastic region ( Inm) followed by a rapid drop in force. The instability is associated with 
bond breaking, but the bonds rapidly reform as the sphere passes through the capsid wall 
and into the interior. This indicates the potential importance of factors that are not included 
in continuum models. Both the yield elastic stiffness and yield stress are about an order 
of magnitude higher than experiment, but vary logarithmically with rate. The yield stress 
often varies logarithmically with rate in thermally activated systems [3], but it is unusual 
to find this variation in the stiffness. This suggests that the simulation rates are higher 
than the structural relaxation rate of the virus and that new behavior may occur at lower 
velocities. 

In this paper we consider a model that is less computationally intensive, but includes 
coarse-grained structural information. The model is used to explore the effect of molecular 
structure on nanoindentation of CCMV and cowpea mosaic virus (CPMV), whose symmetry 



is that of a rhombic tricontahedron 
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23(1 . Capsids in both systems contain of order 300 000 



heavy atoms, which makes it difficult to study large mechanically-induced conformational 
changes with all-atom simulations. Instead we use structure-based coarse-grained molecular 
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dynamics models that take the experimentally determined native structure as an input to 



derive the effective Hamiltonian in a phenomeno 



the detailed description presented in refs. |19|-|21| and this approach has been success: 



ogical way. Our implementation builds on 



ul in 



describing both folding and mechanically- induced unfolding of individual proteins [19|, |21|. 
The coarse-graining reduces the number of effective protein atoms by an order of magnitude. 
Treating solvent molecules implicitly through random force and damping terms leads to even 
greater economies by eliminating atoms and increasing the effective time step. 

We find that the mechanical response of the molecular model of CCMV is qualitatively 
consistent with experiments and continuum calculations, with a long elastic region followed 
by a sudden drop in force. CCMV and CPMV have similar radii, widths and densities and 
thus would be expected to have similar responses in continuum models. In contrast, the 
coarse-grained model predicts that CPMV is an order of magnitude stiffer and breaks at a 
higher force and smaller displacement. These large changes in mechanical response result 
from small changes in bonding, suggesting that the mechanical response of virus capsids 
may be highly variable. Experiments on CPMV would provide a test of this prediction. 

The local deformations in the molecular and continuum models 12| for CCMV are com- 
pared directly. One significant difference is that the molecular model undergoes a gradual 
and reversible rotation in the elastic region. More strain is also accommodated near the 
walls in the molecular model. Dramatic differences between the models are observed at the 
end of the elastic region where the response becomes irreversible. In the continuum model, 
nonlinear it ies are associated with buckling. In the molecular model, nearly all of the inter- 
protein bonds break at the end of the elastic region. This breaking transition appears to be 
thermally activated and its location depends on trajectory. The mean force and displace- 
ment at rupture decrease with decreasing deformation rate, because there is more time for 
thermal activation. 



THE MODEL SYSTEMS 

CCMV is a member of the bromovirus group of the Bromoviridae family. It consists 
of a single-stranded RNA molecule encapsulated by a 180-protein icosahedral capsid. The 
native CCMV capsid is stable around pH 5.0, where it is characterized by the outer and 
inner diameters of 286 and 195 A respectively [2^. At pH 7.0 and low ionic strength, 
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the capsid adopts its swollen form in which its linear sizes are increased by about 10 %. 
Nanoindentation studies at pH 5 yield a spring constant of about 0.15 N/m 5|, |6j. At pH 
6, which is half-way to the situation in which the swollen forrn_^arises, the spring constant 
falls by about a factor of about 3 compared to the native form [6]. Furthermore, there is no 
longer a sudden drop in force that would indicate mechanical instability. 

The symmetries of the native capsid are those of the fullerene Ceo- The proteins form a 
closed shell that has the symmetry of a truncated icosahedron comprised of 32 faces (making 
90 edges and 60 vertices). 12 of these faces are pentagons and 20 are hexagons. A pentagon 
consists of five icosahedral asymetric units ('kites') and a hexagon contains six such units. 
Molecularly, a pentagon corresponds to a pentamer and a hexagon to a hexamer. There are 
thus altogether 180 units, each containing a single molecule referred to as a chain. Even 
though the chains are sequentially identical, one distinguishes three kinds of chains: A, B, 
and C. The pentamers contain only the A-type chains whereas the hexamers have three 
chains of type B and three chains of type C. All chains contain 190 amino acids in a /3-barrel 
fold, but their conformations have not been determined fully through X-ray crystallography. 
In particular, the locations of the first 41 amino acids in chain A and the first 26 in B and 
C are not known. These undetermined segments of the proteins are presumably in more 
random or mobile regions that do not contribute to the scattering or structural integrity. In 
our model, we make use only of the known parts of the structure and the whole capsid is 
represented by 28620 atoms |2J,|25j. 

The conformation of the swollen form of CCMV at pH 7 has been determined through 
cryo-electron microscopy but at a much poorer resolution than the structure at pH 5. 
The resolution is comparable to the separation between atoms that defines the structure 
in our coarse-grained model. We constructed a model based on this structure, but found 
that it was not mechanically stable even at T = 0. This could imply that atoms that 
are misplaced or missing from the structure are crucial to stability. Note that the virus 
disassembles at slightly higher pH and may be nearly unstable at pH 7. We are not aware of 
indentation experiments at pH 7. Given the large reduction in measured stiffness between 
pH 5 and pH 6 [6] and the fact that swelling of the structure mainly occurs above pH 6, 
it would be interesting to extend the experiments to the fully swollen regime. Structural 
information at pH 6 would also allow us to compare our model to experiments in the partially 
swollen state. 
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23j. The 



CPMV belongs to the picornavirus superfamily and can be represented 
triacontahedron with 30 rhombic faces that make 60 edges and 32 vertices 
capsid is assembled from two viral proteins - subunits S and L - which form three /3- 
sandwich domains, A, B, and C in an asymetric unit. They contain 189, 187, and 182 amino 
acids respectively. Domains B and C form hexamers and each L subunit comprises B- and 
C-type domains. The A-domains form pentamers. Altogether, CPMV contains 60 units 
corresponding to the total of 33 480 C° atoms. In contrast to CCMV, the positions of all 
amino acids are determined. 

The angle- averaged radial distributions of the C" atoms in the two capsids in their native 
states are shown in Figure [TJ The average radial position is 119.6 A in CCMV and 124.3 A 
in CPMV. The distributions of radii are also comparable, with total widths of about 55 A 
and rms variations of 11.13 A and 10.54 A, respectively. In continuum theory the stiffness 
associated with compressing a thin shell scales as h"^ / R where h is the thickness and R the 
radius jsj. Given the similarities in the number and radial distribution of C" atoms one 
would expect similar stiffnesses, but we will see that this is not the case. 

The coarse-grained models of the capsids generalize our previous approach for large con- 



formational changes in single proteins 
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3l-|2l|. We take the native coordinates of all heavy 
atoms in all capsid proteins from the structure files and assign van der Waals spheres to the 
atoms (enlarged by a factor that takes into account attraction 26|). If the spheres belonging 
to amino acid i, as counted along the "global" sequence, overlap the spheres related to amino 
acid j then we declare existence of a native contact between amino acids i and j. Most of 
the native contacts are between acids in the same protein, but there are also contacts that 
link separate proteins. Contacts corresponding to the i,i-|-2 interaction within a protein are 

are much weaker 

than the hydrogen bonds. Altogether, there are 62 426 native contacts in the model CCMV 
and 90 420 contacts in the model CPMV. This 45% increase is much greater than the change 
in the number of C" atoms or their radial distribution, and may be an important factor in 
explaining why we find a much greater stiffness for CPMV than CCMV. 

The interactions between amino acids that form native contacts are described by the 
Lennard- Jones (LJ) potential V{r) = 4 e[(^^^ — {jr-^ ], where is the distance between 
the C°'s in amino acids i and j, and aij is determined for each pair to ensure that the 
potential minimum coincides with the experimentally measured native distance. The binding 
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energy parameter e can be determined by comparing theoretical and experimental data for 



(see also ref. 



stretching of single proteins 27|. Our latest estimate 28| is that e is of order 110 pNA 



29|). Interactions along the backbone of each protein are described with a 



harmonic potential with spring constant 50e/A^. This is almost an order of magnitude stiffer 
than the LJ interaction and ensures that the backbone bonds do not break. 

The interactions between atoms that do not form native contacts are represented by a 
purely repulsive potential obtained by truncating the Lennard- Jones potential at its min- 
imum at 4 A. The potential number of non-native contacts is of order 5x10*. To avoid 
calculating all these separations at each step, we make a Verlet list jsol. The list comprises 
all C" atoms that are less than 12 A apart from another atom and is made every 2000 
time steps. 

Langevin noise and damping terms are introduced to provide an implicit solvent 31 |. 



The amplitude of the noise is determined by the temperature, T. All simulations described 
below were done at ksT = 0.3e, where kB is the Boltzmann constant. Estimates of e for our 
model place ksT/e between 0.3 and 0.4 at room temperature. We find similar behavior for 
indentation over this range of temperatures although there are small changes in quantitative 
values. The model is simulated in the overdamped limit and the characteristic time scale r 
is set by the damping rate 2m/ t in the Langevin thermostat, where m is the average mass 



33j. The 



of an amino acid. For proteins in water r is of order 1 ns as argued in refs. [32, 
time step used to integrate the equations of motion with a fifth order predictor-corrector 
algorithm [aoj is r/200. 

Nanoindentation is implemented by placing the capsid between two repulsive parallel 
plates. This geometry corresponds to the limit of a tip with large radius and has been 
considered in continuum calculations. The experimental geometry is often a flat substrate 
with a tip whose radius is of the same order as the virus, but continuum results show that 
this has little effect on the stiffness. The plate potential scales as h~^^ 34], where h is 
the separation between the plate and an amino acid. In the initial state, the plates are 
far enough apart that they do not interact with the capsid. The plates are then brought 
together by increasing the speed of both plates symmetrically to a value Vp over 2000 r. In 



most cases Vp = 0.0025 A/r giving a combined speed of 0.005 A/r that was found to be slow 



enough to produce quasistatic results in studies of protein stretching 
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23] . We will also 



present results at twice this velocity and half this velocity and refer to these as fast and slow 
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compression, respectively. While these velocities (~ 500 /i/s) are higher than experimental 
velocities (0.1 to 10 fim/s), they are slow enough that stress can be transmitted across the 
capsid before the separation has changed substantially. 

We determine total force from the capsid on each plate and then average the magnitude 
of these forces to obtain the total compressive force that would be measured by an AFM. 
The inset in Figure |2] shows that deviations between the forces on the two plates are similar 
to temporal fluctuations in each. Larger deviations with longer duration are observed if the 
velocity is too high for stress equilibration. In this and subsequent figures, the mean force 
is plotted as a function of the plate-to-plate separation, s. While the plate potential has 
a finite range, it changes so rapidly that the effective space available to the capsid is only 
slightly smaller than s. 

Since capsids are not spherical, their elastic properties should depend on their orientation 
relative to the compressive axis. We consider two orientations here. Most results are for the 
case where one of the symmetrically equivalent 2-fold axes isnerpendicular to the plates. 



This corresponds to the z-axis as defined in the structure file [251 ■ The second situation has 
a 3- fold axis normal to the plates. This corresponds to the [1,1)1] direction in the structure 
file and is denoted by [1,1,1]. Note that these orientations represent the initial configuration 
and the capsids can rotate under the applied load. The observed rotations are relatively 
small and are discussed further below. 

RESULTS 

Force-Separation Curves 

Figure [2] shows an example of a nanoindentation trajectory for the CCMV capsid. No 
force is observed for separations greater than 284 A. The force then rises linearly as s 
decreases to about 175 A. This corresponds to a compression by about 35% from the point 
of first contact, a surprisingly large range of Hookean response, but comparable to what 
is observed in experiments on capsids. Note that experimental data are typically plotted 
vs. the displacement of the tip towards the substrate, which becomes more positive as s 
decreases. One must also be careful to remove the compliance of the AFM cantilever in 
comparing to experimental numbers. 
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At the end of the Hookean regime the force drops rapidly by about 40%. A comparable 
drop is observed in experiments . At still smaller separations the force rises rapidly. 
Given the observed width of the radial distributions in Figure [H it is not surprising to find 
that the top and bottom of the virus are in direct contact at these small separations, leading 
to strong repulsive forces. 

In the regime where the force rises linearly, the nanoindentation process is nearly re- 
versible: An inversion of the direction of motion of the plates results in a force which 
returns along the same straight line. Once the capsid has been compressed to or beyond the 
sudden drop in force, the deformation is irreversible. To illustrate this behavior we reversed 
the velocity of the walls after compressing to 50, 100, 125, 150, 175 and 200 A and the results 
are plotted as dotted lines in Fig. |21 The behavior is only reversible for s=200 A. The small 
drop at s just greater than 175 A is already signalling the onset of plastic deformation that 
leads to extremely irreversible behavior. The nature of this plastic deformation is discussed 
further below. 

Figure [3] shows that nanoindentation of a capsid with the [1,1,1] orientation produces very 
similar results. The diameter along this orientation is slightly bigger, so the force starts at 
a larger s, but values of and k are very similar to those in Figure [2l The onset of 
irreversibility also coincides with the end of the Hookean regime. Note that while the capsid 
was free to rotate, the capsid remained in roughly its initial orientation in both Figures [2] 
and [HI 

The velocity dependence of our results is illustrated in Figure HI Three different realiza- 
tions of the system are shown for the velocity considered above, as well as for velocities that 
are slower and faster by a factor of two. The different realizations at a given speed exhibit 
essentially the same behavior at large and small separations, but the onset of irreversibil- 
ity occurs at different forces and separations. This variability suggests that the force drop 
represents an activated transition from the stable linear regime to another state. The rate 
dependence is also consistent with this conclusion. In particular, thermal fluctuations have 
more time to activate the transition as the velocity decreases, and thus the average force at 
the transition decreases with decreasing velocity. The nature of the transition is discussed 
further below. 

The stiffness of the capsid in the linear regime is relatively insensitive to velocity, with 
fluctuations comparable to variations with speed. The results are consistent with a very 
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small (5%) increase in stiffness with speed. A weak linear dependence is expected since 
atoms are able to relax more closely to the free energy minimum as the velocity decreases. 
A stronger velocity dependence is observed at very small s. This is the regime where the 
top and bottom of the capsid are in hard sphere contact. The rapid rise of the repulsive 
force with decreasing separation makes the system more sensitive to unrelaxed positions in 
this regime. 

The CPMV capsid is found to behave in a qualitatively similar way. As shown in Figure |5l 
there is a regime where the force rises linearly, followed by a sudden drop in force. The force- 
separation curve is reversible in the early stages of compression, but becomes irreversible 
after the end of the linear regime. The main differences between the capsids are that CPMV 
is much stiffer in the linear regime and fails at a higher yield stress and smaller degree of 
compression. The peak force Fm ~ 15 e/A is a factor of three larger than for CCMV, while 
the effective spring constant, k ~ 0.5 e/A^, is an order of magnitude higher. Changing 
the orientation of the capsid to [111] produces similar results (Figure E]). The peak force is 
slightly larger, but the difference is within the range of fluctuations between different runs. 

The large difference in stiffness of the two viruses is surprising. Continuum theory predicts 
that k = 2.25Eh^ / R where E is the elastic modulus, h the shell thickness and R the radius 
{5I. As shown in Fig. [H CCMV and CPMV have similar radii R ~ 120 A and thicknesses 
h ^ 20 A. They also have the same interaction energies within our model, and similar 
numbers of amino acids. Yet the measured stiffnesses would imply elastic moduli that differ 
by an order of magnitude: 0.007 eA^^ and 0.07 eA^^ for CCMV and CPMV, respectively. 
To put these values in perspective, they are one to two orders of magnitude smaller than 
the elastic modulus of a Lennard- Jones crystal with the same interactions E=55 e/a^ = 1.2 
e/A^. This large reduction must reflect the fact that the number of bonds constraining the 
motion of each segment of the protein is relatively small. 

Atoms in an fee crystal have z = 12 nearest neighbors or 6 bonds per atom. CCMV 
and CPMV have only 2 or 3 native contacts per atom, respectively. Including the neighbors 
along the backbone of the protein, this corresponds to about 6 neighbors for CCMV and 8 for 
CPMV. The minimum number of neighbors required for stability in 3 dimensions is z = 4, 
and studies of rigidity percolation indicate that the elastic modulus is a strong function z. 
It seems likely that this explains why increasing the number of native contacts per atom by 
only 45% (and total bonds per atom by 30%) can increase the effective stiffness of CPMV 



10 



by an order of magnitude. The observation that small changes in protein binding can have 
a dramatic effect on stiffness suggests that further experimental studies may reveal large 
variations between capsids. We are not aware of any nanoindentation studies of CPMV, 
but it would be interesting to see if the trends observed in our model are reproduced in 
experiment. 



Changes in bonding 

The onset of irreversibility in both CCMV and CPMV is associated with rupture of 
native contacts rather than a buckling transition like that observed in elastic shell models 
I12I . Figure |8] shows the total number, ric, of unbroken native contacts (rjj < 1.5(jjj) 
as a function of s. In the initial stage of compression, the number of unbroken contacts 
remains essentially constant. Thermal fluctuations cause a very small number ( 0.1%) of 
bonds to fluctuate in and out of contact, but the structure remains intact. This behavior 
extends to the end of the Hookean regime where deformation becomes irreversible. At this 
point Tic drops sharply and then saturates as the virus is compressed further. About half 
the contacts break during this transition in both CCMV and CPMV. 

As seen in the force plots described in the previous section, the onset of irreversibility 
depends on thermal fluctuations and tends to occur at larger s for slower compression. The 
drop in Uc always coincides with the onset of irreversibility and shifts in the same way. The 
dotted line in Fig. [8] illustrates an earlier onset of bond rupture in CCMV at a lower rate. 
Note that the values of ric for the different rates coincide before and after rupture occurs. 
Only the point at which thermal fluctuations lead to bond breaking changes. 

Figure [8] shows the total number of contacts, but it is interesting to examine intra and 
inter-protein bonds separately. For the case of CCMV there are three times as many in- 
traprotein contacts (not counting bonds along the backbone of the protein) as interprotein 
contacts: 47166 vs. 15300. When bonds begin to break, the number of broken intrapro- 
tein bonds is slightly higher, but the fraction of broken bonds is dramatically different. By 
s = 149 A over 95% of the interprotein bonds are broken, while more than half the in- 
traprotein bonds remain. This striking result indicates that rupture of the capsid occurs 
mainly along the joints between proteins, something that would be difficult to capture with 
homogeneous continuum models. 
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Interprotein bonds have the same energy parameters as intraprotein bonds, but equihb- 
rium simulations also provide evidence that the joints where proteins meet are weaker. The 
mean squared variation in the length of native contacts in equilibrium should be given by 
the thermal energy divided by an effective spring constant that is influenced by bonds to 
surrounding atoms. For CCMV, the ratio of the mean squared fluctuation to the equilib- 
rium bond length is more than 80% larger for interchain interactions than for intrachain 
interactions, implying that the interprotein regions are roughly half as stiff as intraprotein 
regions. This greater fluctuation indicates that there are fewer reinforcing contacts to limit 
deformation and prevent rupture along the borders of proteins. 

In CCMV, proteins fold into kite shaped pieces and three kites (A, B, and C chains) 



combine to form triangles that tile the surface 



24l |. There are no intrachain bonds connecting 



the kites. In CPMV, three trapezoids combine to form triangles. Two trapezoids on adjacent 
triangles are formed from a single molecule, providing direct intrachain links [2^]. Given 
the mode of rupture, it seems likely that these intrachain bonds between triangles play an 
important role in the high stiffness and rupture force of CPMV. 



Structural Changes in CCMV 

In this section we focus on structural changes in the CCMV capsid during deformation, 
using the trajectory corresponding to Figure [2] to illustrate the changes. Figures l9l fTOl and 
[TT] provide different representations of the conformation of the CCMV capsid at various 
stages of indentation. Figure |9] shows all C° atoms with different gray scales for different 
groups of proteins. Snapshots show the evolution of the capsid shape from equilibrium to 
one of the first states with a significant number of broken bonds, s = 164 A. The walls first 
flatten and compact the region near the plates, with the size of the flattened region growing 
as s decreases. The capsid remains relatively undeformed in the central region between 
walls. 

Figure [10] shows a closeup of the 2000 atoms that are closest to the top plate. In the 
native state, the top of the virus is relatively rough. In the initial stages of deformation 
{s = 264 A) a few of the highest atoms are pushed down but only a small fraction of the 
capsid surface feels the hard repulsive potential from the walls. The force is also very small 
in this range of s, and the stiffness may be slightly smaller than at lower s (Fig. [2]). The 
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number of atoms pushed up against the hard wall potential grows as s decreases. By s = 239 
there is a pronounced flattened region and parts of all proteins shown are in contact with 
the surface. 

Large structural changes such as buckling and rotation are difficult to identify in all-atom 
renderings. As noted above, the capsid is constructed from triangular clusters containing 
three protein chains. A skeleton representation of the structure can be obtained by following 
the center of mass of the three chains forming each triangle. As shown in Figure [111 the 
centers of the triangles lie at the vertices of a truncated icosahedron in the native state. 
In the (110) orientation, a line connecting two vertices is at the top and bottom of the 
virus. Note that this representation makes it appear that there is a sharp edge at the top 
and bottom, which would make the configuration appear unstable against rotations. The 
all-atom rendering shows that the actual surface is fairly spherical and there is no sharp 
edge. 

There are two pentamers and two hexamers that meet at the edges at the top and bottom 
of the capsid. As s decreases from 284 A to 239 A these faces become nearly horizontal, and 
by s = 214 A many of the associated vertices have become coplanar. There is also a slight 
clockwise rotation of the capsid that is discussed below. Further compression to s = 189 A 
seems to produce buckling with the top (bottom) edge pushed below (above) its neighbors. 
These structural results are quite similar to those from finite element calculations with a 

n 

continuum model based on the mass density of the capsid [12]. This model shows that 
the curvature at the top and bottom of the capsid begins to change sign for compressions 
greater than 70 A, which corresponds to s = 214 A. There is a pronounced region of inverted 
curvature at a displacement of 100 A, which is close to the case of s = 184 A in Figure [TTJ 
It is also interesting that for this capsid orientation the continuum calculation finds little 
change in the derivative of force with respect to displacement as buckling occurs, while there 
is a pronounced change for other orientations. Note that we did not observe softening before 
bond rupture in either (100) or (111) orientations. 

For s < 164 A bond breaking becomes important in our simulations. Fig. [TT] shows that 
the capsid has lost most of its structure by s = 139 A. Atoms from the top and bottom of the 
capsid make direct contact at these small separations, leading to a large hard core repulsive 
force on the plates. Bond breaking is not included in the continuum calculations and one 
conseqeunce is that they predict a gradual decrease in the slope of force-displacement curves 
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as s decreases to 164 A, while our calculations show a precipitous drop in force at this point. 

Figure [12] illustrates the spatial distribution of broken bonds for the simulation shown in 
Figure [TTl Each native contact is characterized by its initial "height" h, or more precisely, 
the z-coordinate of its geometrical center in the native conformation. The top panel shows 
the height distribution of native contacts. The bottom panel shows the number with each 
initial height that are broken at various plate separations. Bond breaking initiates just 
before s decreases to 164 A, which is just at the end of the linear portion of the force- 
displacement curve in Figure |2l Breaking initiates near the bottom of the virus for this 
trajectory, but failure initiates at the top in other trajectories. As seen in Fig. [11] and 
continuum calculations, stresses and strains are largest in these extremal regions. 

At s = 164 A about 2000 contacts are broken. Decreasing s by just 5 A ruptures nearly 
all of the bonds that eventually fail in the bottom region and initiates failure at the top of the 
capsid. A further decrease by 5 A has spread the failure over most of the capsid. The number 
of broken bonds then saturates at lower s as shown in Figure [8] The saturated distribution is 
illustrated by the results for s = 139 A in Figure [12] Note that the distribution mirrors the 
number of native contacts, indicating that a fairly constant fraction of the initial bonds has 
broken throughout all regions of the capsid. As noted above, almost all the bonds binding 
different proteins together have failed, leading to complete loss of the capsid structure. 

We now focus on quantifying the distribution of deformation within the capsid, by con- 
sidering the 2;-displacements of individual atoms. They can be illustrated directly through 
a scatter plot in which the atoms' ^-coordinates at separation s, h{s), are plotted against 
their native z-position, h. Figure [13] shows such scatter plots for the C" atoms belonging 
to selected pentamers (3 out of 12) and hexamers (3 out of 20). In each case the n-mers 
are chosen to represent the top, center and bottom and to prevent overlap of points in the 
scatter plot. Results are shown for s = 202 A and s = 164 A where bond breaking becomes 
evident. 

At both plate separations the central n-mers are nearly undeformed. The final height 
differs from the initial height (solid line) in two ways. One is a shift up or down that reflects 
the small clockwise rotation mentioned in regard to Figure [TT] The second is a small change 
in slope. For both plate separations the slope is 0.93 for the pentamer, while for the hexamer 
the slope decreases from 0.96 to 0.93 with decreasing s. The decrease in slope from unity 
results mainly from compression, but is also affected by rotation of the n-mers. The mean 
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rotation of the capsid is only about 4° for s = 202 A. 

The behavior of the outer n-mers is very different. In each case, the largest values of \h\ 
are very close to the height of the walls, s/2. As these outer atoms are pushed inwards, 
they displace the entire thickness of the capsid with them. This leads to a broadening of 
the scatter plots by an amount comparable to the capsid thickness. Note that the envelope 
corresponding to the largest values of \h\ and the outer edge of the capsid does not merge 
smoothly with the lines for the central n-mers. This reflects a rapid change in the surface 
normal at the edge of the central n-mers that is also evident in Figure [111 

To see if the deformations of the capsid in our simulations are consistent with continuum 
models, we compared them to the finite-element results of Gibbons and Klug jl2|. Since the 
finite element nodes do not correspond to individual atoms, horizontal slices of width 1 A 
at selected heights were chosen and nodes at the inner and outer edges of the capsid were 
followed. For the finite-element model these are the inner and outer nodes of the mesh and 
for the molecular model they are the inner and outer C° atoms. The continuum calculations 
moved just one wall so that the center of mass moved. We have corrected for this so that 
zero height corresponds to the center of mass in both cases. 

Figures [H] and [15] compare the finite-element and molecular models for two plate sep- 
arations. Figure [14] shows molecular model results for s = 202 A which corresponds to a 
compression of Z = 82 A. These results are compared to finite element results at Z = 75 
A. Figure [T5] presents results for Z = 116 A in the continuum model and Z = 120 A in the 
molecular model, which corresponds to the onset of bond breaking (s = 164 A). 

For both separations the molecular results show a spread in final heights for atoms with 



the same small initial height 



361]. This reflects the symmetry breaking rotation of the capsid 



that is evident in Figures [TT] and [13] Symmetry breaking is suppressed in the continuum 
model by the lack of thermal fluctuations and the fact that a friction force was imposed 
to suppress sliding of the virus relative to the wall. As shown in Fig. [TS] the spread in 
results for the molecular model results from coherent shifts in the height of each n-mer 
through a relatively rigid rotation. For each n-mer there is a linear relation between initial 
and final height. The continuum model shows a linear relation for all nodes with slightly 
different slopes for the inner and outer atoms. A least mean squares fit to heights of all 
nodes between ib40 A gives 0.90 for the continuum model at Z = 75 A. The same fit for 
atoms in the molecular model gives a slope of 0.94, indicating that there is significantly less 
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strain in the central region. 

The deviations between continuum and molecular models are slightly larger near the top 
and bottom of the capsid. Both show a spread between inner and outer edges, because the 
outer surface is in contact with the wall and the inner surface remains at a fixed distance 
below the outer surface. The continuum model gives a smoother curve for the outer surface, 
while the outer height appears to vary non-monotonically in the molecular model. A local 
dip in final height at an initial height of around 100 A is visible for both wall separations. 
The region where the final height corresponds to contact with the wall also extends over 
a larger range of heights in the molecular model. This is consistent with less strain being 
accommodated in the central region, and more in the regions near the walls. 

The radial strain for both models was also calculated. Defining f as the initial position 
relative to the center of mass and dr as the change in this vector, the radial strain = 
r ■ drjr ■ r. Figure [T6] contrasts the molecular and continuum response for s = 202 A and 
Z = 75 A, respectively. In both cases, the strain is negative at large \z\ where compression 
has fiattened the top and bottom. In the central region the capsid bows out and er is 
positive. For the continuum case the inner and outer regions form smooth parabolas that 
cross at intermediate z. There is some symmetry breaking due to the onset of buckling 
on one side. In the central region, the radius expands by an average of about 15%. The 
molecular model shows substantially less expansion, with an average of less than 5%. The 
envelope of the data is also less smooth, showing a nearly constant expansion in the center 
and a rapid change in slope for > 60 A. This is consistent with the apparent change 
in slope in Figure [I3] and may refiect localized bending at joints between proteins in this 
region. Note that Fig. [15] shows that the fiattened region of the capsid has extended down 
to atoms with initial heights of order 60 A by s = 164 A. It may be that this leads to failure 
because strain can no longer be accommodated by the outer regions. 

To illustrate the development of deformations with decreasing separation, the center of 
mass height z and rms variation Az was calculated for atoms comprising each of the n- 
mers. Figures [T7] and [18] show results for all pentamers and hexamers, respectively. The 
locations of the centers of mass indicate a fairly smooth transition to a fiattened "sandwich" 
state as the wall separation becomes comparable to the thickness of capsid walls. At large 
s, the height of the top and bottom n-mers changes much more rapidly than those in the 
center, indicating yet again that strain is concentrated in the regions near the walls. The 
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mean displacement of n-mers at the same intermediate height shows a relatively smooth 
splitting, that is associated with the symmetry breaking rotation seen before. The splitting 
is not sudden and is reversible when the load is removed. The mean rotation is about 4° for 
s = 202 A. 

Note that one of the n-mers in each pair stays at a nearly constant height while the other 
n-mer moves linearly towards the center. In contrast, the n-mers that start in the center 
move away from the center with symmetry related pairs displacing by equal and opposite 
amounts. Once significant bond-breaking sets in, the heights change rapidly. Note that the 
relative height of outer and more central n-mers changes sign, indicating buckling. It is also 
interesting that the centers of the outer n-mers become equal to the centers of more central 
n-mers at the point where bond breaking sets in. These heights are also equal to the height 
where rapid bending is evident in Fig. [161 

The trends in the rms height change mirror the shifts in the mean. Initially there is a 
fairly linear increase in Az with the greatest slope for the n-mers farthest from the equator. 
The rms displacements are very large, comparable to the mean motion in many cases. This 
reflects the rotation of the outer n-mers as they move into contact with the wall. When bonds 
begin to break, all of the n-mers show rms variations of 20 to 40 A which is comparable to 
the thickness of the capsids and to the plate separation at the smallest s. The n-mers are 
made up of several proteins that are no longer bonded together at small s and thus deform 
incoherently. 

SUMMARY AND CONCLUSIONS 

A coarse-grained molecular model was used to study the mechanical response of CCMV 
and CPMV. Both show a linear elastic response at small deformations followed by a rapid 
drop in force. The response is reversible in the elastic regime and becomes irreversible 
after the force drop due to rupture of bonds within the capsid. Although the two capsids 
have similar radii, widths and densities, their mechanical response is dramatically different. 
CPMV is an order of magnitude stiffer and fails at a much smaller compression and much 
larger peak force. 

The dramatic differences between CCMV and CPMV are the result of a relatively modest 
change in the number of bonds. Insight to this dependence comes from comparing the 
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response to elastic shell models. While the effective modulus from fits to CPMV is an order 
of magnitude larger than that for CCMV, it is more than an order of magnitude smaller 
than that of an fee crystal with the same bonding energies. The crystal would also fail at 
much smaller deformations. The mean number of neighbors drops from 12 in the fee crystal 
to 8 in CPMV and only 6 in CCMV. Four are required to provide mechanical stability in 
3 dimensions and studies of rigidty percolation show that the modulus is a strong function 
of the number and distribution of neighbor bonds. Thus viruses may be able to alter their 
stiffness and ability to accommodate large strains through small variations in the number 
of interprotein bonds. It would be interesting to perform experiments on CPMV and other 
virus capsids to see if the effective modulus does vary by an order of magnitude and whether 
this mechanical property is correlated to the environment or biological functions of the 
different capsids. 

Elastic continuum models produce a nonlinear drop in force when the shell buckles. Since 
there are no thermal fluctuations, buckling occurs at a well-defined displacement where the 
system becomes linearly unstable. In the molecular model the force drop is associated with 
rupture of almost all the interprotein bonds, allowing the capsid to fragment. The capsid 
did not reform on the time scale of our simulations, but experiments show some recovery 
over times of order minutes jsf. 

The simulation results are consistent with the force drop being an activated transition 
where breaking of a few bonds leads to stress transfer and a cascade of additional bond 
breaking. In particular, different trajectories at the same deformation rate follow a common 
force curve before and after the transition, but the transition occurs at different displace- 
ments. As the wall separation is decreased more slowly, there is more time for thermal 
activation and the transition occurs at lower forces and displacements. Experiments show 
a similar trend in the onset of irreversibility with velocity [5], while there is no shift in the 
location of the buckling transition in deterministic continuum models. As is typical of ac- 
tivated transitions, the stiffness in the elastic region is relatively insensitive to deformation 
rate [l8]. An unusual activated rate dependence of the stiffness was observed in recent all- 
atom simulations with a much smaller probe and higher velocities [l7|. As pointed out by 
the authors of this study, the rate dependence they observe suggests that their deformation 
rates were too high to allow the virus structure to relax completely. 

Quantitative comparison to experiments is complicated by rate dependence and uncer- 
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tainties in the binding energy within our molecular model. The latter also leads to uncertain- 
ties in the appropriate temperature. Our slowest rates are two to four orders of magnitude 
faster than experimental rates and the wall displacement that produces irreversibility ranges 
from 150 to 90 A. The range of rates and statistics are too small to extrapolate these re- 
sults to the displacements of about 40 A where irreversibility is observed in experiments. 
Continuum models give buckling at roughly twice this displacement, but buckling does not 
produce a strong drop in the force, particularly for the (100) orientation. The molecular 
model gives buckling at about the same displacement as the continuum model and also shows 
that it does not have a strong impact on the stiffness. The sharp force drop in the molecular 
model occurs at larger deformations when bonds begin to break. At the rates considered 
here, bonds break after the buckling transition, and the onset of breaking shifts relatively 
rapidly with velocity. One might expect that buckling produces highly strained bonds that 
are more susceptible to activated rupture and that the velocity dependence might become 
much slower when the transition occurs at the smaller displacements found in experiments. 

Experiments on CCMV at pH=5 give a peak force of = 0.60 ± 0.04nN and a stiffness 
(-dF/ds) of = 0.15 ± 0.01 N/m in the elastic region [5^. In our simulations, k ~ 0.05 
e/A^ and the peak force is about 4 e/A at the lowest velocity. Using a binding energy of 
e = 110 pNA, these values correspond to A; ~ .055 N/m and ~ 0.4 nN, which are of the 
same order of magnitude as measured values. It is not surprising that agreement is better 
for Fm since e was obtained from fits of simulated force curves during mechanical unfolding 
to experiments that, as here, were done at lower rates. Increasing e would increase both 
Fjji and k, but lowering the rate would bring Fm back toward the experimental value and 
also decrease the displacement at failure toward the 40 A found in experiment. A larger 
e would also decrease the effective temperature, which also increases k and Fm- Thus the 
quantitative agreement is as good as can reasonably expected from a coarse-grained model. 

The local deformations within the CCMV capsid were compared to continuum calcula- 
tions. While the continuum model captures the gross features of the deformation there are 
deviations that must reflect variations in the local elastic properties due to the distribution 
of native contacts. In the molecular model, more of the deformation is accommodated in 
the top and bottom of the capsid. The central region is compressed less, and also shows 
substantially less radial expansion than the continuum model. There is a relatively large 
change in angle at a height of about 60 A that suggests bending. Given the other evidence 



19 



that the regions between proteins are less rigid, it seems hkely that bending is locahzed 
between proteins. More detailed studies of local deformations and stresses is in progress. 
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FIG. 1: The radial distributions of the atoms in the CCMV (dotted line) and CPMV (solid 
line) capsids. The solid squares indicate the mean radii and the horizontal lines show the sizes of 
the standard deviations. 



and all n-mers are considered now. 



FIGURE CAPTIONS 
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FIG. 2: Force on confining plates from CCMV capsid atoms as a function of the plate-plate sep- 
aration in one trajectory corresponding to Wp=0.0025 A/r. The solid line corresponds to nanoin- 
dentation. The dotted lines show the force observed when the plates move back at the same speed 
after reaching the separation indicated by detachment of dotted lines from the solid line. The inset 
shows the separate forces acting on the top (solid) and bottom (dashed) plates. The average of 
these forces is the solid line in the main figure. 
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FIG. 3: Similar to Figure[2]but for nanoindentation along the [1,1,1] direction. 
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FIG. 4: Variation of the force on a CCMV capsid with separation for different trajectories and 
compression speeds. The solid lines in the center correspond to three trajectories at Vp = 0.0025 
A/r. The three dotted and dashed lines correspond to trajectories at half and twice this speed, 
respectively. 
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FIG. 6: Similar to Figure[5]but for nanoindentation along the [1,1,1] direction. 
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FIG. 7: Similar to Figure U] but for the model CPMV. Two trajectories are shown for Vp = 0.0025 
A/r (solid lines) and single trajectories at half (dotted line) and twice (dashed line) this speed. 
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FIG. 8: The number of native contacts that remain unruptured as a function of s for CCMV (soUd 
Hne) and CPMV (dashed hue) in trajectories corresponding to those shown in Figures [2] and [5l 
respectively. A native contact is considered unruptured if the two amino acids are within a distance 
of 1.5 (Tjj. The dotted hne is for compression of CCMV at half the speed and corresponds to the 
lowest trajectory in Figure [H 



30 




FIG. 9: The structure of the CCMV capside with all C locations shown at s=284 A (top left), 239 
A (top right), 189 A (bottom left, and 164 A (bottom right). The smallest value of s corresponds 
to the start of the nonlinear regime. 
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FIG. 10: Snapshots of the 2000 highest C" atoms for s equal to 284 A, 264 A, and 239 A, from 
top to bottom. All of these snapshots are in the Hookean regime. 
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FIG. 11: The vertex representation of the CCMV capsid at various stages of indentation charac- 
terized by the indicated values of s. The sohd circles are closer to the viewer and the open circles 
are in back. 
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FIG. 12: The top panel shows the z-dependent distribution of the native contacts, nc{h), as 
explained in the text. The bottom panel shows the distributions of contacts that are broken at the 
indicated values of s in A. 
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FIG. 13: Scatter plots of final heights of C" atoms as a function of initial heights in the native 
structure. The top two panels are for s=202 A and the bottom two panels are for s=164 A. 
Hexamers and pentamers were selected to represent the polar and equatorial regions. 
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FIG. 14: Scatter plots of final vertical positions as a function of initial height for outer (solid 
symbols) and inner (open symbols) boundaries of the capsid. The top panel shows postions of 
atoms for s=202 A, and the bottom panel shows no de p ositions from the continuum model of 
Gibbons and Klug at a wall displacement of Z = 75 A. [12]. 
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FIG. 16: Scatter plots for the radial strain as a function of height corresponding to a wall separation 
of s = 202 A. The top panel is for the molecular model and the bottom from the data provided by 
Gibbons and Klug from a continuum model [l^. 
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FIG. 17: The bottom panel shows the height of the centers of mass of atoms belonging to 12 
individual pentamers. The top panel shows rms variation in the height change of atoms belonging 
to the same pentamer. The outer lines in the bottom panel correspond to the polar pentamers and 
the lines in the middle to the equatorial regions. The same symbols are used in the top panel. 
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FIG. 18: Similar to Figure [T71 but for the 20 hexamers. 
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